/*
 *  simKernel.cu
 *  sparticlus
 *
 *  Created by Lita Gratrix on 4/16/10.
 *  Copyright 2010 __MyCompanyName__. All rights reserved.
 *
 */
 
#include "simKernel.h"

#define kNumberOfArraysPerCalc 3

//Helper function that checks for a cuda error after a cuda call. If an error
//has occured, print the error message and exit sparticlus
void checkCUDAError(const char *msg){
    cudaError_t err = cudaGetLastError();
    if( cudaSuccess != err) 
    {
        fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString( err) );
        exit(-1);
    }                         
}

//CUDA kernel that calculates the new paricle position (posPrime) based on the initial position (posArray) and velocity (velArray)
__global__ void calculateNewPosition(int fps, float incr, float *posArray_d, float *velArray_d, size_t pitchInFloats, unsigned int partCount, float *posPrimeArray_d){

	int idx = blockIdx.x * blockDim.x + threadIdx.x; //id the row of the executing thread
	int idy = blockIdx.y * blockDim.y + threadIdx.y; //id the column of the executing thread
	
	//if the thread index is greater than the array index, don't bother processing the rest of the kernel 
	if  (idy >= kVArraySize || idx >= partCount){
		return;
	}
	
	posPrimeArray_d[idx * pitchInFloats + idy] = posArray_d[idx * pitchInFloats + idy] + ((velArray_d[idx * pitchInFloats + idy] / fps) / incr);
	
	//Debug print statements
	/*
	printf("posArray_d[%d][%d] = [%f],", idx, idy, posArray_d[idx * pitchInDoubles + idy]);
	printf(" velArray_d[%d][%d] = [%f], ", idy, idx, velArray_d[idx * pitchInDoubles + idy]);
	printf(" posPrimeArray_d[%d][%d] = [%f]\n", idy, idx, posPrimeArray_d[idx * pitchInDoubles + idy]);
	*/
}

//set up the cuda structures required to calculate new position data on GPU
void simKernel(int fps, float incr, float *posArray, float *velArray, unsigned int partCount, float *posPrimeArray_h){
	int i;
	//Debug print statements
	/*printf("Upon entering simKernel, incr is %f, and partCount is %d.\n", incr, partCount);
	 
	 for (i = 0; i < partCount * kVArraySize; i+=3){
	 printf("posArray[%d] is: [%f,%f,%f]\n", i/3, posArray[i], posArray[i + 1], posArray[i + 2]);
	 printf("velArray[%d] is: [%f,%f,%f]\n", i/3, velArray[i + 0], velArray[i + 1], velArray[i + 2]);
	 }
	 */
	
	//pointers to arrays on the GPU
	float *posArray_d, *velArray_d, *posPrimeArray_d;
	
	//pointer for pitch (width in bytes that cudaMalloc assigns per row) of 2D arrays. 
	// I am using the same pointer for multiple allocations becuase they all have the 
	// same size so should all return the same pitch value; some risk here
	size_t pitch = 0; 
	
	int memSize = sizeof(float) * kVArraySize; //size of width in bytes required by a 3 element array of floats
	
	//allocate memory on the GPU for 2D array of vectors with partCount number of rows
	cudaMallocPitch((void **) &posArray_d, &pitch, memSize, partCount);
	checkCUDAError("allocating device memory for posArray_d");
	cudaMallocPitch((void **) &velArray_d, &pitch, memSize, partCount);
	checkCUDAError("allocating device memory for velArray_d");
	cudaMallocPitch((void **) &posPrimeArray_d, &pitch, memSize, partCount);
	checkCUDAError("allocating device memory for posPrimeArray-d");
	
	//copy array data from host to device
	cudaMemcpy2D(posArray_d, pitch, posArray, memSize, memSize, partCount, cudaMemcpyHostToDevice);
	cudaMemcpy2D(velArray_d, pitch, velArray, memSize, memSize, partCount, cudaMemcpyHostToDevice);
	checkCUDAError("copying array data to device");
	
	int pitchInFloats = pitch / sizeof(float);
	
	//Debug print statements
	/*
	 for(i = 0; i < partCount; i++){
	 printf("posArray_d[%d] = [%f, %f, %f].\n", i, posArray_d[i * pitchInFloats + 0], posArray_d[i * pitchInFloats + 1], posArray_d[i * pitchInFloats + 2]);
	 }
	 printf("sizeof(float): %d\n", sizeof(float));
	 printf("sizeof(char): %d\n", sizeof(char));
	 printf("memSize: %d\n", memSize);
	 printf("pitch: %d\n", pitch);
	 */
	
	//dimension the kernel attributes
	int nBlocks = ceil((float)partCount/ (float)kBlockHeight); //number of blocks required 
	dim3 threadsPerBlock(kBlockHeight, kVArraySize); //nubmer of threads per block - 2D array
	
	//kernel call
	calculateNewPosition <<< nBlocks, threadsPerBlock >>>(fps,  incr, posArray_d, velArray_d, pitchInFloats, partCount, posPrimeArray_d);
	cudaThreadSynchronize();
	checkCUDAError("Kernal execution");
	
	//copy new position from device to host
	cudaMemcpy2D(posPrimeArray_h, memSize, posPrimeArray_d, pitch, memSize, partCount, cudaMemcpyDeviceToHost);
	checkCUDAError("cudaMemCpy");
	
	//Debug print statements
	/*
	 for (i = 0; i < partCount * kVArraySize; i+=kVArraySize){
		printf("posPrimeArray_h[%d] is: [%f, %f, %f]\n", i/kVArraySize, posPrimeArray_h[i], posPrimeArray_h[i + 1], posPrimeArray_h[i + 2]);
	 }
	*/
	
	//free device memory
	cudaFree((void *)posArray_d);
	cudaFree((void *)velArray_d);
	cudaFree((void *)posPrimeArray_d);
	checkCUDAError("freeing device memory");

	return;
}

